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Abstract. We develop perturbation theory and physically motivated resummations 
of the perturbation theory for the problem of a tracer particle diffusing in a random 
media. The random media contains point scatterers of density p uniformly distributed 
through out the material. The tracer is a Langevin particle subjected to the quenched 
random force generated by the scatterers. Via our perturbative analysis we determine 
when the random potential can be approximated by a Gaussian random potential. We 
also develop a self-similar renormalisation group approach based on thinning out the 
scatterers, this scheme is similar to that used with success for diffusion in Gaussian 
random potentials and agrees with known exact results. To assess the accuracy of this 
approximation scheme its predictions are confronted with results obtained by numerical 
simulation. 



PACS numbers: 05.20.-y, 66.10.Cb, 66.30.Xj 
1. Introduction 

The bulk transport properties of quenched random media play an important role 
in engineering sciences such a hydrology JLj. Calculational techniques to evaluate 
bulk transport properties such as effective diffusivities, mobilities, conductivities and 
permeabilities thus have far reaching applications. In systems such as zeolites, disorder 
modifies the diffusion of reacting chemical species and thus plays a crucial role in 
diffusion limited chemical reactions |2]. The study of such systems is also relevant 
to the evaluation of short-time transport properties in non-quenched random media 
where the time scale on which a tracer diffuses over the typical correlation scale of the 
advecting field is much shorter than the time scale controlling the temporal evolution of 
the advecting field. Systems with quenched disorder may also be used as paradigms for 
systems with temporal disorder and may sometimes be used to obtain bounds on the 
transport properties of non-quenched fields. 

Recently there has been much interest in the nature of the glass transition, in 
particular as to whether it is a purely dynamical or a static transition. Consider a 
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tracer particle diffusing and interacting with other identical particles as could be the 
case in a liquid or gas. If inertial effects are not important we can assume that the 
motion of the particle is controlled by over-damped Langevin dynamics. In the liquid, or 
ergodic, phase we expect the tracer particle to have a non-zero self diffusion constant n. 
Now consider a situation in which the same particle diffuses in a quenched background 
where all the other particles have been frozen in a particular configuration, possibly 
chosen from Gibbs-Boltzmann equilibrium ensemble. Let the diffusion constant in this 
quenched system be Kg. It has been shown [3^ that Kg < /t; a result that seems physically 
plausible since, if the background particles can move about, the cages which could trap 
the tracer in the quenched background, will break up on some time-scale and free it 
to disperse more quickly than in the quenched case. This consideration recently lead 
Szamel |1] to study this problem as a test of the mode coupling theory for interacting 
particle systems. 

A Langevin particle diffusing in a random potential is often invoked to describe, at 
a suitably coarse grained level, the transport properties of certain random media. Often 
the first choice of statistics for the potential is Gaussian. If the potential has a finite 
correlation length, then we can subdivide the system into uniform blocks whose size is 
larger than that correlation length. Consequently the diffusion processes within each 
of these blocks will be independent of one another. The particle will stay within each 
block for a random time r, that will be distributed similarly but independently for each 
block. For Gaussian statistics we expect the mean value of r to be finite. Thus, viewed 
as coarse grained hopping process between blocks, the motion of the tracer will appear 
to be an uncorrelated random walk and thus asymptotically exhibit normal diffusion. 

The calculation of the diffusion constant of a Langevin particle in a Gaussian 
potential, specified therefore solely by its two point correlation function A(x), has been 
attempted via a wide variety of approximation schemes 0. The most successful to 
date is the self-similar renormalisation group (RG) scheme of 0. Other approximation 
schemes such as self-consistent perturbation theory, and the Hartree-Fock approximation 
perform poorly when confronted with known exact results and numerical simulations. 
The modification of the diffusion constant in the presence of inertial effects |^ can 
also be well explained using a Gaussian closure scheme, within the Martin-Siggia-Rose 
path integral formalism [2] . This closure approximation is rather reminiscent of certain 
aspects of mode-coupling theory |H] and is related to the dynamics of mean field spin 
glasses [Oj (the approximation becomes exact in the limit of infinite spatial dimension 
and when the field strength is suitably scaled). This Gaussian closure predicts a glass 
transition [Hj (signalled by a vanishing of the diffusion constant and ageing in correlation 
and response functions) for a Langevin process diffusing in a short range correlated, finite 
variance Gaussian random potential, however from our previous discussion we do not 
expect such a transition in finite dimensions. 

In this paper we shall consider the diffusion constant of a Langevin tracer particle 
diffusing in a system of quenched random scatterers which are independently and 
uniformly dispersed through out the volume. This corresponds to the situation found 
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in zeolites where the scatterers are usually charged ionic sites. We draw attention to 
the fact that this model is distinct from that of the Lorentz gas in which an undamped 
Newtonian particle is scattered by hard spheres in three dimensions or hard discs in two 
dimensions (overlapping or non-overlapping) [TU] . 

We develop a double expansion for the self diffusion constant in interaction strength 
and the density of scatterers and determine the range of validity of the Gaussian 
approximation for the random field statistics. We then evaluate the density expansion 
to first order and use renormalisation group ideas to resum this expansion. This strategy 
leads to a result which agrees with known exact results for diffusion the presence 
of scatterers in one and two dimensions We identify the situations where this 

approximation appears to be accurate by comparison with numerical simulations of 
diffusion in the presence of hard and soft sphere scatterers. A distinct trend appears 
from our study, a renormalization group calculation based on the weak density expansion 
appears to work very well in problems where the scatterers attract the tracer particle, 
but is considerably less performant where the scatterers are repulsive. 

A number of useful mathematical results for Langevin particles diffusing in random 
potentials, which we shall make use of in this paper, are given straightforward derivations 
in the appendix, thus making the paper self contained. 



2. The model 



The position Xj, of the Langevin tracer particle subject to a white noise, ri{t) and a 
force generated by a potential 0(x), satisfies the stochastic differential equation 



dt 

where 



v^r7(t) - AV0(Xt) , (1) 



{7],it)7],it')) = 26,,5it-t') , (2) 

and (■ ■ ■) denotes an average over the white noise. The Einstein relation implies that 
the local, or bare, diffusivity k and the coupling to the potential gradient A are related 
by the equation 

where T is the absolute temperature in appropriate units. The probability density, 
p(x, t), for the position of the particle obeys 

^ = nV^p + AV ■ (pV0) (4) 

= - Hp , (5) 

where equation (0) defines the positive operator H. The effective diffusivity, Kg, 
describing the dispersion of the particle at large times and distances is defined by 

n, = — \im^ (6) 



Perturbation theory for the effective diffusion constant in a medium of random scatterersA 



where D is the dimension of space. The mean squared displacement of the particle is 
calculated via 

(X,2) = |dxxVx,t) . (7) 
In this paper we will consider potentials of the type 

N 

0(x) = ^r(x-x„) , (8) 

n=l 

generated by scatterers frozen in the volume V. Here x„ is the position of scatterer 
n and in what follows all the x„ will be taken to be independently and uniformly 
distributed in the volume V. This situation where the scatterers are uniformly and 
independently distributed could arise in circumstances where the medium is formed by 
the rapid quench of a very high temperature hquid or gas to a solid phase, the scatterers 
then becoming frozen in the high temperature configuration. The density of scatterers 
is p = N/V . The thermodynamic limit is V ^ oo at fixed p . 

The potential V{x.) for a scattering particle at the origin will be taken to be 
rotationally invariant that is ^(x) = l^(|x|). The quenched disorder in this problem is 
therefore completely specified by the positions, x„, of the scatterers and the disorder 
average is given by 

{0), = l[^ fd^^ O . (9) 

The Green's function for the problem is G'(x, x') = i/^^(x, x'). It satisfies the 
equation 

kV^G' + AV ■ (G"V0) = -5(x - x') . (10) 
The Fourier transform G'(k, may be defined to be 

G'(k) = j dy exp(-«k ■ Y)G\y + x', x') . (11) 



Because the system is statistically homogeneous in space we may for most purposes set 
x' = in the above without loss of generality. From equation ()10|) we have in Fourier 
space 

'5w = ;s^-;s5/(lF''"'^('"«'(k-q, . (12) 

The effective diffusion constant Kg can now be obtained from the long distance behaviour 
of the Green's function [.5. . In Fourier space we have 

G'(k) \-r as |k| ^ . (13) 

Kgk^ 

One can calculate Kg exactly for this problem in one dimension [TTl and one finds 

Ke = Kexp (^-2p J dx (cosh(/5y(x)) - 1)) . (14) 

Remarkably the expression ()14|) is independent of the sign of jS. For positive (3 the 
particle is attracted to the minima of the potential (p and for negative /? it is attracted 
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Figure 1. Diagrammatic expansion for G"(k) 



to the maxima. To disperse however, the particle must visit both minima and maxima, 
due to the absence of saddles in one dimension, and so one could argue that, in some 
sense, both minima and maxima were equally important for determining the dispersion. 
The result is still however rather intriguing from a physical point of view. Another 
interesting property of the one dimensional result is that if one considers a system with 
two types of scatters, types a and b, one type of density pa with potential Va and the 
other of density with potential VJ, one finds that 

Ke{a)Ke{b) 

Keia + b) = , (15) 

where Ke{a + b) is the effective diffusivity for system with both scatterer types and Ke(a) 
and Ke{b) are the effective diffusivities for the system with only particles of type a and 
b respectively. 

An exact result was also recently found in two- dimensions. If one considers a system 
of charged scatterers where 

N 

4>{x) = 1nV{^ - X„) , (16) 

n=l 

and where g„ is the charge of the scatterer at Xn which takes the (quenched) value 
g„ = ±1 with probability 1/2 so that the system is statistically electro-neutral, then it 
was shown that 

Kg = /texp ^— p J (ix [cosh (/9V(x)) — 1]^ . (17) 

The interesting feature in this case is again that the multiplicative structure seen in 
equation (jlSp in the case of two families of independent charged scatterers still holds. 



3. Weak disorder expansion 

In the non-neutral case in two and higher dimensions one must revert to approximate 
methods to calculate He- We consider the perturbation expansion generated by iterating 
equation (fT^ and which is shown diagrammatically in figure (^. 

The Feynman rules for this perturbation theory (before averaging over the disorder) 

are: 
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• The total momentum at each vertex is conserved. 

• Each sohd horizontal line with momentum k flowing through it is the bare 
propagator l/nk^. 

• Each wavy line with momentum q flowing into it carries a factor — Ak- q</)(q) where 
k is the incoming (from the horizontal left) momentum. 

• Each wavy line momentum q in the diagram is integrated with the measure 
dq/(27r)^. 

One now has to average over the disorder to obtain 

Gik) = {G'ik)), , (18) 

where k = |k|, as the averaged Green's function is clearly isotropic. We now note that 
the Fourier transform of is given by 

0(q) = E^(q)exp(-^q-x„) . (19) 

n 

Therefore each vertex contains n terms of the form ^(q) exp(— zq ■ x„). If one considers 
a diagram coming from the expansion of all the sums at each vertex containing the 
particle index i, k times at the vertices z/i • • ■ z/^, then the contribution depending on the 
position Xj of scatterer i is given by the prefactor 

k 

P = exp(-z5]x,-q,„) . (20) 

m=l 

Now, the disorder average with respect to x^ yields 

{p)^=^-^5ij:^^j ■ (21) 

m=l 

Hence the total momentum flowing up into vertices with the same scatterer index is 
zero. We also note that in the unaveraged Feynman rules any vertex with upward 
flowing momentum q,^ = causes the diagram containing it to be zero. In each diagram 
one has N choices for the scatterer index at the flrst vertex. Every distinct scatterer 
in each diagram carries a factor of 1/V from averaging over its position. A diagram 
with r distinct scatterer indices at its vertices (but with the positions of each index 
flxed) therefore carries a factor because there is a factor for the number of ways 
to choose the scatterers plus a factor of 1/V from the averaging over the r distinct 
scatterer positions. We can represent this ordering graphically as is shown in flgure 
(j21) where the one, two and three vertex contributions are shown, the flrst column 
represents the flrst scatterer chosen in the flrst vertex, the momentum at each blob is 
the momentum flowing into that vertex. We continue this development where the u^'^ 
row corresponds to the vertex. In each of these diagrams the sum of the momentum 
down each column is zero due to equation ()21|) . Also, from the Feynman rules, any 
diagram with a blob containing zero momentum must be zero; along with momentum 
conservation this means that any non-zero diagram must contain each scatterer zero 
times or at least twice. The diagram flgure (j2t) is thus zero as are all but the flrst of 
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Figure 2. Blob diagrams to order 
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Figure 3. Non-zero blob diagrams of order A* 

the diagrams in figures ^Bp) and (I2fc)- The order of each diagram written this way is 
p'^ X X^, where c is the number of columns, corresponding to distinct particles, and r is 
the number of rows, corresponding to the number of vertices. 

In figure we show the non-zero four-vertex contributions: figure(jnk) has only one 
column (scatterer) and figures ^jp), (jSfc) and represent the two scatterer diagrams. 

We note that for large p the dominant diagrams with an even number 2n of rows 
are those where each column has only two blobs and that this diagram, in units where 
K = 1, is of order (pA^)". Note that these diagrams are formed by pairing the vertices 
and reproduce the diagrams for a particle diffusing in a Gaussian random field [H] , where 
the field ijj has the correlation function given by 



(^(q)V'(q')). = (27r)V'5(q + q')l^(q)y(-q) 



(22) 



For diagrams with an odd number of vertices 2n + 1 the dominant diagrams in p 
are those with n — 1 pairs of blobs with a single column containing three blobs and the 
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order of these diagrams is, in units where k = 1, p"-i_\2n+i _ {^pX^^'^Xj p_ Thus in the 
hmit A — > whilst keeping pA^ = constant, one recovers the hmit where a system of 
scatterers generates a potential which can be treated as Gaussian for the purposes of 
calculating Kg- This limit can also be seen in exactly soluble one dimensional and two 
dimensional models 

The modified Feynman rules leading to systematic expansion of G(k) in terms of A 
are thus: 

• At order A*^ draw all distinct blob diagrams with r rows. 

• Each blob diagram carries a factor of p'^ where c is the number of columns. 

• Each blob diagram has a corresponding Feynman diagram with the same rules as 
before except it is V^(q) at each vertex instead of 0(q). 

• The total momentum down each line is zero and each independent remaining 
momentum is integrated with the measure c?q/(27r)-^. 

Finally the series can be Dyson resummed in terms of one-particle irreducible 
diagrams by writing 

where S is given by the summed one particle irreducible blob diagrams multiplied by 
K^fc^, that is to say the corresponding diagram with the first in going and the last out 
going bare propagators amputated. As usual a one-particle irreducible diagram is one 
that does not factorise. All the diagrams of order p are clearly one-particle irreducible. 
In figure ©, figure (jSh) is clearly reducible and the others are irreducible. 
In terms of the one particle irreducible diagrams 

Ke = K — lim —TIT- ■ (24) 

In what follows we shall use the fact that V^(q) = V{q) where q = |q|, due to the 
radial symmetry the the potential V . Using the perturbation theory rules developed 
here we find that to order (3'^ 

pp^ f dq p(3^ f dq dp ( (p-q)2\ ~ ~ ~ 

7 = ^- ^y(2^^^^)+:^y (2^ (2^ V - -p^ 1 ^(^)^(^)^(ip + ^1) 



2 /q4 r _ 7 n 2 



2^2 



(27r) 



D 



+ p^pR . (25) 



The principal difficulty with the above perturbation theory is the identification of the 
O(fc^) term in each one-particle irreducible diagram. The term R not written explicitly 
in the above suffers from this problem and it is not clear how the resulting integral 
should be expressed for efficient evaluation. Also expansion as a power series in f3 is 
presumably not very accurate where Kg deviates significantly from its bare value, as the 
first term is 0(/?^). 
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4. Weak density expansion 

One of the principal technical difficulties with the weak disorder expansion as presented 
in the previous section is that its convergence can be expected to be poor, and it may 
possess divergences for singular potentials. In addition, as the interested reader may 
verify, the identification of He by extracting the terms 0{k'^) in the expansion of S(A;) 
is algebraically rather difficult. In the case of scatterers it is natural to examine the 
expansion of Kg in terms of the density p of the scatterers. The small density expansion 
was pioneered by Maxwell [T^ in the context of effective dielectric constants. We shall 
develop this analysis for the problem of the scatterers in a compact and direct manner. 
We start with the formal expression for the diffusion constant which is given in terms 
of the Green's function by 

^ ~ ~ Uz 

where Z is the one particle partition function for a particle in the potential given by 



j d^d^ G"(x, x) V(/)(x) ■ V(/)(x') exp (-/3(/)(x')) , (26) 



Z= / rfxexp(-/30(x)) . (27) 



V 

This result does not appear to be widely known (or at least used) in the physics 
literature. It can be derived from static considerations ^2] or via the stochastic 
differential equation for the diffusion process For a proof based on the stochastic 
approach we refer the reader to the appendix of this paper. 

The expression equation ()26|1 may be used to reformulate the weak disorder 
expansion by exploiting the perturbation expansion for G' . We see that to order p the 
effective diffusion constant is the result of independent contributions from scatterers, 
that is 

— = I-EttI- / G'i(x,x',x„)V0i,„(x)V0i,„(x')exp(-/50i,„(x')) , (28) 

where t6p denotes Kg calculated to ffist order in p. In the above (/)i^„(x) = y(x — x„) is 
the potential due to the scatterer at x„ and Zi „ the corresponding one particle partition 
function. We have denoted by Gi(x, x',x„) the diffusion Green's function for a tracer 
particle in the presence of only one scatterer at x„. Note that without loss of generality 
we may take all the x„ = and in the limit of large V and iV, when l^(x) goes to zero 
sufficiently quickly as |x| ^ oo so that the integral 



F = y rfx' [1 - exp (-/3V(x'))] (29) 

is finite, then 



-p- 



. ^ /c/xrfx Gi(x,x',0)VV(x)-Vr(x)exp(-/3V(x')) .(30) 

K D J 

We note here that to order p the value of Kg is independent of the nature of the 
correlations between the x„ and the formula equation ()30|) relies only on the fact that 
the distribution of scatterers is homogeneous and isotropic. We now define 



w x 



(3k J dx.'Gi{x,x.',0)VV{x.')exp{-(3V{x.')) , (31) 
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and by the definition of Gi, the Green's function for diffusion in tlie presence of a single 
scatterer taken to be at tlie origin, we liave tliat 

VVi(x) + /3V ■ Wi(x)Vy(x) = -PViVi^) exp (-/3\^(x)) . (32) 

Using tlie spherical symmetry of the problem we write w = exp{—f3V{r))g{r):s./r where 
r = |x| and g obeys 

rf2 D-l d D-1 ^d d ^d 

—g + —g g-p g V + P—V = . 33 

dr^ r dr dr dr dr 

In terms of g we thus have 



|-exp(-/?V(r)) 



9ir) , (34) 



K D 

where Sd is the area of a unit sphere in D dimensions. The behaviour of g at large r 
needs to be determined. Note that if V decays sufficiently quickly at large r then we 
can write equation (jH^ as 

VV(x) ^-/5ViV^(x)exp(-/3r(x)) , (35) 

for large r. The effect of the source term must die away as r ^ cxd for rapidly decaying 
potentials and so at large r we may write w = V'?/' where ip obeys 

= exp (-/3l^) - 1 , (36) 

which has the solution 

V^(x) = j dx'AB^(x, x') [1 - exp (-/5r(x'))] , (37) 

where 

is the inverse of the Laplacian in D dimensions (the case D = 2 is to be interpreted as 
a logarithm). Now at large |x| 

w(x) ^ V^(x) ^ C^-^ j dx' [1 - exp (-/5K(x'))] (39) 



when the integral F defined in equation ()29|) is finite. For F to be finite we must have 
that V{r) decays quicker than to ensure convergence as r ^ oo. In this case we 
see that equation ()39|) then implies that g{r) ~ l/r^~^ for large r. 

The boundary conditions on g are that 5f(r) as r ^ oo (for D > 2) from the 
previous discussion and, if the potential V is sufficiently well behaved, we expect that 
g{0) is finite. In greater than one dimension equation (j33|l is difficult to solve explicitly. 
However for the case of soft sphere potentials of the form Vs{r) = e for r < a and 
Vs{r) = for r > a, one can solve equation (jH^ explicitly. In this case it is useful to 
rewrite equation (jHSI) as 

(exp{-/3Vir))r''-'-^g]+expi~(3Vir))r'''%l-D)g-r'''^-^exp{-(3Vir)))=0 .(40) 

(Ji I \ (Ji I / (Ji I 
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^('•-'^«)+'-'(l-Z»«^0 . (41) 

for all r ^ a, supplemented by the conditions for matching the discontinuity at r = a 
which are 

gia-) = g{a+) (42) 

a^~V^(a+)-exp(-/3e)^(a-))-a^-Vl-exp(-/3e)) = . (43) 
dr dr 

The general solution to equation ()41|) is 

^W = ^ + ^'^ ' (44) 

and the boundary conditions at r = and as r — > oo along with the continuity and 
jump condition across r = a give 

1 - exp(-/5e) 

q[r) = -r ;r < a 

' D - 1 + exp(-/5e) 



l-exp(-/3e) /.rN 
^ ^ a (-)^ ^ ;r > a. (45) 



D — 1 + exp(— /5e) r' 



We thus obtain 



^ = l + ^a^-i(l-exp(-/5e))^(a) 

^ _ pSpa^ (1 -exp(-/3e))^ 

D (D-l+exp(-/3e)) ■ ^ ' 

In the limit of hard spheres, that is to say as e — oo, we obtain 

^ = 1 _ ^^^""^ (47) 
K D{D-l) ■ ^ ' 

This result is divergent at all p for D = 1 as it is clear that in the presence of hard 
scatterers the tracer particle will be confined between its nearest right and left neighbour 
for all time and thus will not diffuse. 



5. Self-similar renormalisation group schemes 

The basic idea of self-similar renormalisation group schemes is to resum the original 
perturbation expansions available to us in a physically motivated way. Consider the 
general problem of diffusion in a random potential where the potential is parameterised 
by A such a way that when A = the potential vanishes and the resulting diffusion is 
free. For example we could use a Fourier based parameterisation where we write 

M^)=f, 7T^0(q)exp(zq-x) , (48) 

Now we define a perturbation 50a (x) to the potential so that 

0A+<5A = 0A + S(j)A ■ (49) 
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In the case of the Fourier parameterisation we have therefore 




(50) 



which is the usual shcing technique in Fourier space. 

One then integrates out the high momentum component 50a perturbatively to find 
an effective theory where the effects of 50a have been taken into account on the running 
diffusion constant denoted by k(A) and the couphng to the remaining gradient field 
denoted by A(A). This is a self-similarity ansatz which means only keeping interactions 
that were in the original problem and hence the only parameters which change are k 
and A (as well as the potential which just has the term 50a removed). One therefore 
has a running diffusion constant k(A) and a running coupling to the gradient field A (A). 
The effective diffusion constant is then given by Kg = k{0). In addition it can be shown 
by general arguments J2] that 



that is to say that the Einstein relation, or fluctuation dissipation relation, is 
satisfied by the renormalised theory at each step of the renormalisation. This 
renormalisation is only valid for the low-momentum component of the remaining drift 
and a possible improvement to the calculation here would be to functionally renormalise 
the remaining field 0a, which amounts to introducing new interactions generated 
by the renormalisation procedure (this approach has been applied to diffusion in an 
incompressible quenched Gaussian velocity field fl]). If V0 = u were a constant applied 
field then A would be the local conductivity coupling the particle to this external field. 

We remark that one may verify from the results on the weak disorder perturbation 
theory in section 2 that the Fourier space based renormalisation group method does not 
reproduce the exact results for scatterering generated potentials known in one and two 
dimensions. The Fourier based renormalisation group does however reproduce known 
exact results in one and two dimensions for Gaussian potentials |3] . However the Fourier 
space slicing of the potential is not only the only decompostion possible. It was noted 
that in the case of diffusion in a Gaussian field a so called t-slicing renormalisation group 
procedure (TH] gave identical results to the Fourier space renormalisation group. In the 
case of anisotropic potentials these two renormalisation group schemes give different but 
numerically very close results. The basic idea behind the t-slicing approach is to write 



where 0' and 0" are two independent Gaussian fields with the same correlation function 
as the original field 0. In the above prescription A is integrated down from one to zero, by 
perturbatively calculating, and averaging over, the effect of term 0". The corresponding 
flow equations for k(A) and A (A) are then integrated. The physical argument for the 
Fourier based renormalisation group scheme is that at length scales 1/A the effective 
diffusion constant and effective 'conductivity' A is determined by fluctuations in the 
field on length scales shorter than this. In the case of the scatterers one can apply a self 



0^ = (A-5A)20' + 5A20' 



(52) 
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similar renormalisation group approach based on thinning out the density and which has 
the advantage of agreeing with known exact results. Concretely, one assume that the 
effect of a density p of scatterers is to generate an effective diffusion process described 
by an effective diffusion constant ^(p) and coupling A(p) on length scales bigger that the 
inter-particle separation. We now throw in more scatterers of the same type of small 
density 6 p. As 5p is small we may use the results of the previous section along with the 
self-similarity ansatz here to write 



6k{p) 6pS 



D 



which integrates to give 

k(p) = Kexp 



dr 



|:exp(-Wr)) 



9[r) 



' pS 



D 



D 



dr 



d 

— exp 
or 



9[r) 



(53) 



(54) 



We see that the effect of the renormalisation group procedure is to exponentiate the 
small p expansion equation (jMj) . Immediately we see from its multiplicative form that 
the result equation ()34|) agrees with the known results in one and two dimensions referred 
to previously. Note it is important that the scatterers are uniformly and independently 
distributed as the added scatterers at each step have no information about those already 
integrated out. Let us also note that on adding an additional Sp of scatterers, the 
characteristic length between these scatterers is /* ~ l/Sp^ which is strictly infinite as 
6p —>■ 0, hence a diffusing particle will, on average, have acquired the effective diffusion 
constant k{p) and conductivity A(p) before encountering a newly added scatterer. For 
a system of soft spheres we have in the RG approach 



Kg 



exD ^'^^^ (l-exp(-/36))- 
^[ D (D-l + exp(-/?e))J • 
Despite its agreement with known exact results in one and two dimensions the 
expression equation (j54|) will not be exact in general. Indeed, if one considers a 
system of hard spheres one expects that there is a percolation transition at a density 
Pc above which the diffusion constant vanishes, the formula equation ()55p is clearly 
insensitive to any such transition; we see that in the hard sphere limit e — oo that 
Ke = exp (^—pSdci^ /{D{D — ^))^- It is inevitable that this approximate result will 
be inaccurate at high densities. In the following section we shall compare the results 
obtained by the various techniques used here with numerical simulations. 

A related problem to the one studied here is the problem of calculating the effective 
diffusivity of a medium with locally varying diffusivity. The diffusion equation for this 
system is 



dt 



[/€(x)Vp) 



(56) 



The effective diffusion constant in this model also corresponds to the effective dielectric 
permeability of a medium with local permeablility k(x) and thus we shall denote it by 
K,^f\ In the case where 



K X 



ftexp(— /30(x)) 



(57) 
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one can relate the diffusion constant for a particle diffusing in the potential at inverse 
temperature f3 with the bare diffusion constant k, to that of the effective diffusion 
constant of the process defined by equation by the remarkably simple formula 

.(^) = ^ . (58) 

We denote by k the spatial average of the local diffusivity or permeability 

K = — /" (ix k(x) . (59) 



VJv 

This result was shown in ^2] via an algebraic approach but it has a simple probabilistic 
derivation which we include in the Appendix. In the case considered here, of uniform 
scatterers, we have 



K = K exp 



P 



j dx(exp(-/5y) - 1) 



(60) 



The study of the effective dielectric constants of a mixture of dielectric spheres suspended 
in an otherwise uniform medium dates back to Maxwell [1^]. If the spheres have 
dielectric constant n' and the background medium has dielectric constant then the 
local dielectric constant may be written as 

/€(x) = K exp (— 0(x)) , (61) 

where 

</) = 5:i^(x-x,) , (62) 

i 

with V{r) given by the soft sphere potential: V{r) = ln(/t/K') for r < a and V{r) = 
for r > a, with a the radius of the spheres. This formulation clearly works when the 
positions of the spheres are chosen to be non-overlapping. However, to order p the 
spheres do not overlap and thus from the calculations developed here, using equation 
PUj) with e = Iii^k/k'), equation (j3Hj) and equation (plj) to order p, we find 

4" pSpa" {i-if . /. , pSoa" fc' \ 

[D — 1) + a 

where c = pS^a^ /D is the volume fraction of the spheres and a = k'/k. In three 
dimensions this result gives 

3c(a-l) , 
^ ~ 1 + ^ 1. 64 

which agrees with Maxwell's result, which he gives as JH] 

^ ^ a + 2 + 2c{a-l) 
K a + 2-c{a-l) ' ^ ' 

to first order in c. Note that the Maxwell result is only accurate to order c and 
equation is a physically motivated resummation of the order c result, it is sometimes 
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Table 1 


. Calculations of nf^ / k 


for 


a — GO 


c 


numerical 


O(c^) result 


Maxwell 


RG 


0.00 


1.000 


1.000 


1.000 




1.000 


0.10 


1.352 ± 0.030 


1.345 


1.333 




1.349 


0.20 


1.821 ± 0.081 


1.780 


1.750 




1.822 


0.30 


2.529 ± 0.168 


2.306 


2.236 




2.459 


U.4U 


o.oyU ± U.ZoU 


2.923 


3.000 




o.ozU 


^0 


4 Qfi7 + 33Q 


3.628 


4.000 
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Table 2 


. Calculations of Hif^ / k 


for 


a = 10 


c 


numerical 


0{c^) result 


Maxwell 


RG 


0.00 


1.000 


1.000 


1.000 




1.000 


0.10 


1.247± 0.011 


1.240 


1.243 




1.252 


0.20 


1.545± 0.028 


1.511 


1.529 




1.568 


0.30 


1.944± 0.046 


1.812 


1.871 




1.964 


0.40 


2.443± 0.046 


2.144 


2.286 




2.460 


0.50 


3.080± 0.053 


2.507 


2.800 




3.080 



referred to in the literature as the Clausius-Mosotti equation. The correction at order 
was calculated many years later by Jeffrey 17]. Despite the fact that the calculation 
carried out in this paper is for spheres which can overlap, it is interesting to compare 
the full renormalisation group result for overlapping spheres with results obtained for 
the effective dielectric constant for a system of non-overlapping spheres via numerical 
calculation. The renormalisation group calculation for overlapping spheres gives 

«:(^') = «:exp(3c^) . (66) 

Brady and Bonnecaze ^H] numerically computed the effective dielectric constant for 
non-overlapping dielectric spheres in a background medium. In Tables (1-3) we compare 
their results with equation (j66|l . the Maxwell formula equation (j65|) and Jeffrey's result 
which is correct to O(c^). For a = oo and a = 10, the RG result is surprisingly close 
to the numerically obtained values and performs better than the Maxwell formula and 
the 0{(?) result. For a = 0.01 it performs less well than the two other results. 

The effective potential for the scatterers in the potential diffusion problem was 
e = — ln(a). Hence if a > 1 then the potential e is negative and thus attractive to a 
particle diffusing in that potential. The diffusion of the particle is clearly slowed down 
by trapping inside the scatterers, and it appears that the RG is picking up this dominant 
effect, correlations between scatterers give a small effect. 
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Table 3. Calculations of nf^/K for 


a ^ 0.01 


c 


numerical 


Uyc ) result 


Maxwell 


RG 


0.00 


1.000 


1.000 


1.000 


1.000 


0.10 


0.859± 0.004 


0.858 


0.859 


0.862 


0.20 


0.729± 0.005 


0.727 


0.731 


0.744 


0.30 


0.600± 0.007 


0.608 


0.614 


0.641 


0.40 


0.504± 0.010 


0.501 


0.506 


0.554 


0.50 


0.406± 0.012 


0.455 


0.407 


0.478 



If we consider the case where a < 1, then the potential e is positive and repels the 
tracer. The slowing down of the diffusion is now mainly induced by another, much more 
subtle, effect: a memory effect which is much more sensitive to the correlations between 
the scatterers. 

6. Simulation Results 

The values of the effective diffusion constant as predicted by Equation and (f^ have 
been compared to the effective diffusion constant measured in Monte Carlo simulations 
of two dimensional systems of soft spheres. The system is a square of linear size L = 20, 
with the spheres of radius, as defined by Equation (jlUj) . a = 1 uniformly distributed in 
the area. As usual periodic boundary conditions are imposed. The simulations consisted 
of 2000 independent runs of a single tracer particle during a total time t^ax = 10^, with 
an elementary time increment 6t = 0.0025. Each particle was evolved in an independent 
realisation of the disorder. The value of was evaluated by fitting the ensemble average 
of X^/t to 4/te + Aln{t)/t at late times (see Equation ()A.12|) ). The errors were then 
estimated from the typical late time fluctuation of the ensemble averaged curve about 
the fitting curve. In Figure (jH) Kg is plotted as a function of T, for p = 0.35 and e = — 1, 
i.e. a case where the scatterers attract, and tend to trap, the tracer. The agreement 
with the RG calculation is extremely good, the RG calculation correctly predicts the 
behaviour of He across a range where it varies almost by a factor of 5. In this case the 
accuracy of the RG is comparable to that when it is applied to the problem of diffusion 
in a Gaussian random potential in three dimensions 

In Figure © i^e is plotted as a function of T, for p = 0.55 and e = 1, thus the 
scatterers repel the tracer. In this case the RG works notably less well and the numerical 
curves appear to be bounded above by the RG and bounded below by the order p result. 

Finally in Figure © Kg is plotted as a function of p for e = oo, the case of hard 
sphere scatterers. Here the perturbative order p result performs extremely well and 
is superior to the RG result. This may be expected as the ratio Xq/hq in this case is 
effectively infinite and the renormalisation in the zero temperature case may not make 
sense. 
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Figure 4. Effective diffusion constant as a function of the temperature, for e = — 1 
and p = 0.35: from numerical simulations (a), low density approximation (b) and 
density RG (c). 



The comparisons here are in accordance with those seen in the last section, the RG 
seems to perform well when the slowing down of the diffusion is due to trapping centres. 
It performs significantly less well in the case of repulsive scatterers where the slowing 
down of the diffusion is due to more subtle correlation effects and presumably needs a 
higher order treatment of terms in the expansion in p. 

7. Conclusion 

In this paper we have considered the problem of a Langevin process diffusing in a 
quenched random potential. The random potential is generated by uniformly and 
independently distributed scatterers. We developed a weak coupling or high temperature 
perturbation theory which can be used to determine the validity of the Gaussian 
approximation to the random field distribution. As expected from the central limit 
theorem, one of the conditions needed to validate the Gaussian approximation is that 
there is a high density of scatterers so that the potential at each point is the sum of 
a large number random variables. However this condition is only sufficient at high 
temperatures. Physically this can be understood as follows: even though, due to the 
central limit theorem, at a randomly chosen point the potential will be Gaussian, at low 
temperatures the diffusion will be dominantly near the minima of the potential which 
by definition belong to the tails of the distribution where the central limit theorem does 
not apply as they are at the extreme points of the distribution. We have also seen that 
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Figure 6. Effective diffusion constant as a function of the density for e = oo 
(overlapping liard spheres): from numerical simulations (a), low density approximation 
(b) and density RG (c). 
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the self-similar RG resummation of the weak disorder expansion does not reproduce 
known exact results. In addition the weak disorder expansion is destined to fail for 
singular potentials, where the potential contains a divergence. We thus formulated the 
self-similar RG resummation in terms of a decimation or thinning procedure in terms of 
the particle density. Direct and indirect comparison with numerical results for hard and 
soft sphere potentials show that the density based RG resummation gives reasonable 
predictions for the bulk diffusivity and agrees with known exact results. It appears to 
work best in cases where the scatterers have a potential which is attractive to the tracer 
and thus act as trapping centres. In the case where the scatterers are not trapping but 
are repelling, diffusion is again slowed down but by a more subtle correlation mechanism, 
this effect is less well described by the RG result. 
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Appendix A. Expression for 

Here we derive the representation equation for the effective diffusion constant. The 
basic idea of the proof in much more mathematical terms may be found in the book 
of Spohn Consider the process in equilibrium in a large volume V at the time 

t = 0. Let the initial position be Xq distributed according to the equilibrium measure 
Peg(xo) = exp(— /50(xo))/2' wherc Z is the one particle partition function normalising 
the distribution. Integrating the stochastic differential equation (H)) between and t we 
obtain 

Xt - Xo = ^/^Bt -xfds V0(X,) , (A.l) 

JO 

where Bj is a standard D dimensional Brownian motion with (B^) = 2Dt. Squaring 
the above and taking its average yields yields 

2DKt = {(Xt-Xoy)+2X{{Xt-Xo)- f' dsV(j){X,))+2X^ f' ds T ds' V0(X,)-V0(X,O .)(A.2) 

Jo JO Jo 

Defining a time dependent diffusion constant of K*{t) = ((X^ — Xq)^) /2Dt, so that 

limt__+oo i^*{t) = f^e this yields 

W = - ^ ((X* - Xo) ■ f ds V0(X,)) -^{f ds r ds' V0(X,) ■ V</.(X,0) .(A.3) 
Ut Jo Ut Jo Jo 

Using the property of detailed balance, which ensures time translation invariance at 

equilibrium, we have for any two functions A and B that 

{A{X,)B{X,)) = (A(X,_,)5(Xo)) , (A.4) 

and the symmetry relation 

{A{Xt)B{X,)) = {B(Xt)A{X,)) . (A.5) 
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Using this it is easy to verify that the second term on the right hand side of equation 
()A.3|) is zero. The third term may be simphfied using 

lit) = {f ds r ds' V0(X,) ■ V0(Xy)) 

JO "'0 

= {f ds r du V0(X„) ■ V0(Xo)) 

JO JO 

= / c/xfixo / ds I c/m V0(x)p(x,xo,m) ■ V0(xo)peg(xo) , (A. 6) 

J JO JO 

where p(x, xo,m) is the transition density to go from xq, to x in time u and we have 
used the fact that the density for the starting point xq is the equihbrium density. We 
now use the fact that formally 

p(x,Xo,m) = exp(-Mif(x,Xo)) , (A. 7) 

to obtain 

/(t) = y"rfxdxo^ dsV(p{^)H-^{l - exp(-s/7)) • V0(xo)peg(xo) 
= t j (ix(ixoV0(x)if~^ ■ V0(xo)peq(xo) 

- j d^d^oV(j){^)H-\l - exp{-tH)) ■ V0(xo)peg(xo) , (A.8) 

where in the above the notation of multiplication by H^^ corresponds to operator 
composition. Using the definition of G' and putting this altogether we find, keeping 
the dominant terms in the large t expansion, that 

r 

Ke = K - — j (iX(iXoV0(x) ■ G'(x,Xo)V0(xo)peg(xo) , (A.9) 

which is the desired result. Note that all the above derivation depends on the fact that 
if is a positive definite operator. In a finite volume this is clearly not the case as there 
is an equilibrium distribution. The result should be understood for long times but such 
that the system has not encountered the boundary of V. If in the above the limit t — > oo 
is taken before |V| — >■ oo one finds that Kg = which is to be expected for a finite system. 
A more detailed discussion of this subtlety can be found in [T^ . 

Another interesting consequence of equation ()A.8|) is that using the zeroth order 
(or free diffusion) result for and exp(— tiJ) in the finite time correction in equation 
()A.8|) . we obtain the expected finite time correction jT^] to the long time diffusive 
behaviour, i.e. 

{X^t) = 2Ket + 0{Vi) -,0 = 1 (A.IO) 
= 4M + 0(ln(t)) ;D = 2 (A.ll) 
= 2DKet + 0{l) ]D>?, . (A.12) 

Appendix B. Relationship between and K^f> 

The diffusion equation equation (j^UI) corresponds to a process which in the Ito 
prescription of the stochastic calculus obeys 



c/X = JKexp(-/?0(X))r7(t)rft - /t/?exp(-/?0(X)) V0(X)rft (B.l) 
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where fi;(x) = Kexp(— /30(x)). We bear in mind that the equihbrium distribution of this 
process can readily be seen from equation (j^Bj) to be flat, i.e Peqi'^) = Now define 
a random time variable 

ds exp{-f3(j){Xs)) (B.2) 



which is clearly a monotonic proper, though random, change of time scales. Clearly we 
may write equation (jB.lj) as 



dX = ^Kexp{f3(f)(X))r]{t)dT - AV0(X)rfr , (B.3) 

and we have used k(3 = X. Formally we may write 

7]{t) = cr/VIt , (B.4) 

where a is a normal random variable of zero mean and variance 2. We can therefore 
write equation ()B.3|1 entirely in terms of the time variable r 

dX = ^r]{T)dT - \V(t){X)dT . (B.5) 

This equation is just that for a particle diffusing in the potential at inverse temperature 
(3. Thus by definition 

(X^) = 2k^Dt (B.6) 

but also 

(X?) = 2K(P)Dt . (B.7) 
Comparing the two expressions which must be the same one finds 

tv^^ —— Kip 



t 



= K,l f ds exp(-/30(X,)) . (B.8) 
t Jo 

We now use the fact that X has a fiat equilibrium distribution to give that (for large 
times) 

i f'ds exp(-/30(X,)) ^ 1 / dxexp(-/30(x)) , (B.9) 
t Jo V Jv 

which thus leads to the desired result equation (j3Hj) . 
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